module TriangleElementManipulation
    use MatrixOperation
    use GaussianIntegration
    use VectorOperations
    implicit none

    integer, parameter :: nNodes = 3
    integer, parameter :: probDim = 2 !< Problem Dimension | Plane Problem

    integer, parameter :: MatrixSize = nNodes * (probDim + 1) !< Heat Fluxes (each dimension) and Temperature problem
	integer, parameter :: XAxis = 1
	integer, parameter :: YAxis = 2

	integer, parameter :: HeatFluxX_Initial = 1 !< Keeping track of Degrees of Freedom index meaning
	integer, parameter :: HeatFluxX_Final   = HeatFluxX_Initial + nNodes - 1
	integer, parameter :: HeatFluxY_Initial = HeatFluxX_Final + 1
	integer, parameter :: HeatFluxY_Final = HeatFluxY_Initial + nNodes - 1
	integer, parameter :: Temperature_Initial = HeatFluxY_Final + 1
	integer, parameter :: Temperature_Final = Temperature_Initial + nNodes -1

	integer, parameter :: Dof_Values = 2
	integer, parameter :: Dof_Try = 1 !< Try values for the element DoF
	integer, parameter :: Dof_Commit = 2 !< Commit values for the element DoF

    !> \brief Triangular Element
    !! \warn The element is not aware of its Neighbours. Is up to the grid decide which is neighbour of which
    !! The Solution of the Discontinuous Finite Element formulation is by Sucessive Substitution, which is equivalent to
    !! iteratively solve a block-diagonal global matrix. This kind of solution is well-suited to parallel environments
    !! with fast interprocess communication systems. Due to this choice of solution procedures the local values of the
    !! Degrees of Freedom must be separated into two categories, the trying value and the committed value.
    type TriangleElement
!// \todo (ronaldo#1#): Modificar Matrix, Dof, OuterVector, NeighbourReNumbering para allocatable e ajustá-los em SetNodesTriangularElement
        real(kind=8), dimension(:,:), allocatable :: Matrix
        integer, dimension(nNodes)                        :: Nodes
        real(kind=8), dimension(:,:)    :: DoF !< Degres-Of-Freedom are comprised for the Try and Commit Values
        real(kind=8), dimension(:,:)           :: OuterVector !< Outer vector oposed to the vertex
        integer, dimension(nNodes, nNodes)                :: NeighbourReNumbering
        real(kind=8) 									  :: S !< Auxiliar Space variable used to make calculations faster
    contains
        procedure, public :: SolutionErrorRange
        procedure, public :: SetNodesTriangularElement
        procedure, public :: SetMatrixTriangleElement
        procedure, public :: SetVectorTriangleElement
        procedure, public :: SetAllNeighboursNumbering

    end type

contains


    !> \brief Returns the Error range estimation for the current Element
    !!
    !! \param this - the element
    !! \return the maximum error over the components of the Degrees-Of-Freedom vector
    !!
	pure function SolutionErrorRange( this ) result (v)
		implicit none
		type(TriangleElement), intent(in) :: this
		real(kind=8) :: v
		v = norm_infty(this%DoF(:,Dof_Commit)-this%DoF(:,Dof_Try))
	end function

    !> \brief Set Nodes Indexes of a Triangular Element
    !! These are the nodes of the grid
    !! \param this - the element itself
    !! \param Nodes - the Indexes of the Nodes on the Grid
    !!
    subroutine SetNodesTriangularElement( this, Node1, Node2, Node3)
        implicit none
        type(TriangleElement), intent(out) :: this
        integer, intent(in) :: Node1
        integer, intent(in) :: Node2
        integer, intent(in) :: Node3


		integer :: statReturn

        this%Nodes(1) = Nodes1
        this%Nodes(2) = Nodes2
        this%Nodes(3) = Nodes3

		allocate( this%Matrix(MatrixSize, MatrixSize+1), stat=statReturn )
		if (statReturn /= 0) then
			print *, 'Unable to allocate memory for Matrix on Triangle Element.'
			return
		end if

        this%Matrix = 0.0d0

        allocate( this%DoF(MatrixSize, Dof_Values), stat=statReturn )
        if (statReturn /= 0) then
        	print *,  'Unable to allocate memory for Degrees of Freedom on Triangle Element.'
        	return
        end if

        this%DoF = 0.0d0

		allocate( this%NeighbourReNumbering(nNodes, nNodes), stat=statReturn )
		if (statReturn /= 0) then
			print *, 'Unable to allocate memory for Neighbour Renumbering Matrix on Triangle Element.'
			return
		end if

		this%NeighbourReNumbering = 0

        allocate( this%OuterVector(nNodes, probDim), stat=statReturn )
        if (statReturn /= 0) then
        	print *, 'Unable to allocate memory for Outer Vector on Triangle Element.'
        	return
        end if

    end subroutine



    !> \brief Groups the individual constructions of the several components of the Linear Systems Matrix
    !!
    !! \param this - The Triangular Element itself
    !! \param
	!!
	subroutine SetMatrixTriangleElement( this, X, ThermCond, Stabilization )
        implicit none
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(nNodes, probDim), intent(in) :: X
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface
		real(kind=8), intent(in) :: Stabilization

		this % Matrix(1:MatrixSize,1:MatrixSize) = 0.0d0

        call SetMatricesE( this%Matrix(1:MatrixSize,1:MatrixSize), X )
        call SetMatricesH_l( this, this%Matrix(1:MatrixSize,1:MatrixSize), X )
        call SetMatricesJ_l( this, this%Matrix(1:MatrixSize,1:MatrixSize), X )

		call SetMatricesE_lr_i( this, this%Matrix(1:MatrixSize,1:MatrixSize), X, ThermCond, Stabilizaiton )
		call SetMatricesH_l_i( this, this%Matrix(1:MatrixSize,1:MatrixSize), X, ThermCond, Stabilizaiton )
		call SetMatricesJ_l_i( this, this%Matrix(1:MatrixSize,1:MatrixSize), X, Stabilizaiton )
		call SetMatrixG_T_i( this%Matrix(1:MatrixSize,1:MatrixSize), X, Stabilizaiton )
    end subroutine


	subroutine SetVectorTriangleElement( this, X, Neighbours, ThermCond, Stabilization )
        implicit none
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(nNodes, probDim), intent(in) :: X
        type(TriangleElement), dimension(nNodes), intent(in) :: Neighbours
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface
		real(kind=8), intent(in) :: Stabilization

		call SetMatricesE_lr_B_i(this, this%Matrix(:,MatrixSize+1), Neighbours, ThermCond, Stabilization )
		call SetMatricesH_l_B_i(this, this%Matrix(:,MatrixSize+1), Neighbours, Stabilization, ThermCond )
		call SetMatricesJ_l_B_i(this, this%Matrix(:,MatrixSize+1), Neighbours, Stabilization )
	end subroutine

	subroutine SetAllNeighboursNumbering(this, Neighbours)
		implicit none
		type(TriangularElement), intent(inout) :: this
		type(TriangularElement), dimension(nNodes), intent(in) :: Neighbours

		integer, intent(in) :: i

		do i = 1, nNodes
			call SetNeighbourNumbering(this, i, Neighbours(i))
		end do
	end subroutine


	subroutine SetNeighbourNumbering(this, n, Neighbour)
		implicit none
		type(TriangularElement), intent(inout) :: this
		integer, intent(in) :: n
		type(TriangularElement), intent(in) :: Neighbour

        integer :: k, f
        do k = 1, nNodes
            f = findNode(k)
            if (f == 0) then
                this%NeighbourReNumbering(n,k) = k
            else
                this%NeighbourReNumbering(n,findNode(this, k, Neighbour)) = k
            end if
        end do
    end subroutine

    pure function findNode(this, n, Neighbour) return (NbIdx)
        implicit none
		type(TriangularElement), intent(inout) :: this
        integer, intent(in) :: n
		type(TriangularElement), intent(in) :: Neighbour

        integer :: NbIdx
        integer :: k
        do k = 1, nNodes
            if (Neighbour%nodes(k) == this%Nodes(n)) then
                NbIdx = k
                return
            end if
        end do
        NbIdx = 0
    end function


	pure function VectorRenumber( this, i ) result (Indexes)
		implicit none
		type(
		integer, intent(in) :: i
		integer, dimension(MatrixSize) :: Indexes

		integer :: k
		do k = 1, nNodes
			Indexes(HeatFluxX_Initial+k-1) = HeatFluxX_Initial+this%NeighbourReNumbering(i, k)-1
			Indexes(HeatFluxY_Initial+k-1) = HeatFluxY_Initial+this%NeighbourReNumbering(i, k)-1
			Indexes(Temperature_Initial+k-1) = Temperature_Initial+this%NeighbourReNumbering(i, k)-1
		end do
	end function

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Calculations Zone Delimiter!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


	!> \brief Calculate the rotation of nodes over a Triangle
	!!
	!! \param n - the edge oposed to vertex n which is going to be calculated the normal vector
	!! \return v - 3-d integer vector with the local indexes of the Nodes
	!!
	pure function rotation(n) result (v)
		implicit none
		integer, intent(in) :: n
		integer, dimension(3) :: v
		select case(n)
			case(1)
				v(1) = 1
				v(2) = 2
				v(3) = 3
			case(2)
				v(1) = 2
				v(2) = 3
				v(3) = 1
			case(3)
				v(1) = 3
				v(2) = 1
				v(3) = 2
			case default
				v(1) = 1
				v(2) = 2
				v(3) = 3
		end select
	end function


    !> \brief Calculate the outer normal vector oposed to de vertex i
    !!	the calculations are preformed using a extension of the problem dimension to generate a orthogonal vector by cross vector product.
    !! \param this - Triangular Element
    !! \param Coordinates - The Coordinates Position of the Nodes
    !!
    subroutine CalculateOuterVector(this, Coordinates)
        implicit none
        type(TriangleElement), intent(out) :: this
        real(kind=8), dimension(nNodes,probDim) :: Coordinates

		real(kind=8), dimension(2) :: Xm
		real(kind=8), dimension(3) :: V1, V2

        integer :: n
		integer, dimension(3) :: p

        do n=1,3
			p = rotation(n)

			!< Xvirt = Xm(extended); Xm(3) = 1.0d0
			!< The projection of Xvirt on the plane is already over the the edge defined by Xp(2) and Xp(3)
			Xm = (Coordinates(p(2),:) + Coordinates(p(3),:))/2.0d0
			V1(1:2) = Xm-Coordinates(p(2),:) !< V1 = XVirt-X2
			V1(3) = 1.0d0

			V2(1:2) = Coordinates(p(3),:)-Coordinates(p(2),:) !< V2 = X3-X2
			V2(3) = 0.0d0

			!< This is already a normal vector
			V1 = Cross(V1, V2) !< Cross vector product is not assumed to have the vectors V1 and V2 well positioned

			Xm = Xm - Coordinates(p(1),:) !< Outer vector of edge oposed to vertex 1
			if (dot_product(Xm,V1(1:2))>0) then
				this%OuterVector(p(1),:) = V1(1:2)
			else
				this%OuterVector(p(1),:) = -V1(1:2)
			end if
			call Normalize(this%OuterVector(p(1),:))
        end do
    end subroutine





end module
